Multipole strength function of deformed superfiuid nuclei made easy 
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We present an efficient method for calculating strength functions using the finite amplitude 
method (FAM) for deformed superfiuid heavy nuclei within the framework of the nuclear density 
functional theory. We demonstrate that FAM reproduces strength functions obtained with the fully 
self-consistent quasi-particle random-phase approximation (QRPA) at a fraction of computational 
cost. As a demonstration, we compute the isoscalar and isovector monopole strength for strongly 
deformed configurations in ^''"Pu by considering huge quasi-particle QRPA spaces. Our approach 
to FAM, based on Broyden's iterative procedure, opens the possibility for large-scale calculations of 
strength distributions in well-bound and weakly bound nuclei across the nuclear landscape. 
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Introduction.- One of the major challenges in the 
many-body problem is the microscopic description of the 
collective motion involving hundreds of strongly interact- 
ing particles. Here, of particular interest is the response 
of the system to a time-dependent external field. In the 
nuclear case, in the small-amplitude limit of nearly har- 
monic oscillations about equilibrium, the phenomena of 
interest include the variety of vibrational modes, and 
characteristic distribution of electromagnetic, particle, 
and beta-decay strength [l|, 0] . 

Most nuclei exhibit strong nucleonic pairing that pro- 
foundly impacts their collective motion 0, |j]. When 
moving away from the stability line, another factor af- 
fecting nuclear correlations, and dynamics is the presence 
of a low-lying particle continuum which provides a vast 
configuration space for nucleonic excitations. Therefore, 
to understand the variety of nuclear modes throughout 
the nuclear chart, a consistent treatment of many-body 
correlations and continuum is required Q . 

This study is devoted to the multipole strength in 
superfiuid deformed heavy nuclei. The theoretical 
method is the quasi-particle random-phase approxima- 
tion (QRPA) applied to the self-consistent configuration 
obtained by means of the nuclear density functional the- 
ory (DFT) Q . QRPA represents a small amplitude limit 
to the time-dependent superfiuid DFT method. In the 
absence of pairing, QRPA reduces to the usual Random 
Phase Approximation (RPA) built atop the Hartree-Fock 
(HF) equilibrium. 

In the electronic DFT, the RPA contribution to elec- 
tron correlation energies has emerged |7|, |8| as an impor- 
tant building block of accurate density functional treat- 
ments of molecules and solids as it combines a number 
of attractive features: it includes the long-range disper- 
sion [3 as opposed to semi-local functionals; it is non- 
pcrturbative and can be applied to small or zero gap 
problems, such as metals [l3| or dissociating H2 111; it 



is nearly exact in the high-density or low-coupling limit, 
and it is parameter-free; finally it is intimately related to 
the microscopic Coupled Cluster doubles theory |12h14| . 
In the nuclear DFT, the applicab ility of (Q)RPA to cor- 
relation energies is more limited |15l4l7| as many nuclei 
have transitional character, i.e., they are close to the crit- 
ical point for the symmetry-breaking where the second- 
order expansion in density fluctuations breaks down [18| . 
In spite of its drawbacks, because of its deep connection 
to DFT, QRPA is still the tool of choice when it comes to 
either spherical or well-deformed nuclei. In addition to 
numerous applications to small-amplitude collective mo- 
tion, QRPA for deformed nuclei may be utilized in the 
calculation of the collective mass for the large-amplitude 
dynamics |19j . 

The advantage of QRPA-I-DFT is that it properly 
takes into account self-consistent couplings giving rise to 
the variety of symmetry-breaking phenomena and, when 
done properly, includes the effects due to the continuum 
coupling. Due to the complexity of the problem, the 
QRPA framework being capable of a fully self-consistent 
description of non-spherical systems has eluded us until 
very recently [l^-Hl. 

A major obstacle preventing the widespread use of 
QRPA has been its high computational cost. In chem- 
istry, this factor has limited applications of this method 
considerably Q- In nuclear physics, most fully self- 
consistent QRPA applications have been limited to spher- 
ical nuclei (see, e.g., [25| and references therein). In spite 
of advanced computational resources available, it is only 
very recently that deformed QRPA calculations atop the 
Hartrce-Fock-Bogoliubov (HFB) equilibrium have been 
carried out H^-il- 

The primary challenge in the conventional matrix 
formulation of (Q)RPA is computation and storage of 
huge matrices of the residual interaction. The recent 
breakthrough came from the realization that both bot- 



tlenccks can be avoided by taking advantage of self- 
consistent DFT solutions and directly employing the lin- 
ear response theory to them (see literature quoted in 
Refs. [IJ, [2^, l27[ ) • Indeed, in both the finite-amplitude 
method (FAM) [2g| and the iterative non-Hermitian 
Arnoldi diagonalization technique l27|, the (Q)RPA ma- 
trix equations are recast into the set of self-consistent 
equations with respect to (Q)RPA amplitudes, which sig- 
nificantly reduces computational effort. The FAM has 
been applied in the RPA variant to Skyrme energy den- 
sity functionals (EDFs) to study giant dipole resonances 
and low- lying pygmy dipole modes [28J, |29| . Recently, in 
its QRPA extension, FAM was used to study monopole 
resonances in a spherical drip-line nucleus ^'^^Sn 30[. 
The iterative Arnoldi diagonalization was first used in 
the RPA variant to electromagnetic strength functions 
in ^"^^Sn [23], and the spherical QRPA extension has re- 
cently been reported [31|. 

While based on the same principle, the actual numer- 
ical implementations of FAM and iterative Arnoldi di- 
agonalization differ. In the applications of FAM, the 
(Q)RPA amplitudes are iterated at desired energies. The 
Arnoldi algorithm is a Krylov subspace iterative method 
for extracting a partial eigenspectrum, i.e., it computes a 
discrete set of states that approximate true eigenvectors. 
(For another nuclear application, an iterative Lanczos- 
based power iteration algorithm for solving the RPA 
equations, see Ref. |32|.') 

The current implementation of FAM has so far been 
done in the coordinate-space representation that requires 
a large number of iterations (in some cases more than 
500-1,000 [2a, [23) to obtain self-consistent amplitudes. 
Here, wc propose a fast and efficient method for solv- 
ing the FAM- QRPA equations in the harmonic oscillator 
(HO) basis using the Broyden iterative scheme [3j, |34| 
previously adopted to HFB equations of nuclear DFT 
[35t . We study the performance of the method and com- 
pare it with the standard QRPA diagonalization method. 
We demonstrate that FAM-QRPA solutions can be ob- 
tained with no more than 40 iterations at modest mem- 
ory requirements of about half a gigabyte for large model 
spaces corresponding to extreme cases of fission isomers 
in the actinidcs. 

Finite Amplitude Method.- The basic formulation of 
the FAM-RPA is presented in Ref. [Hi, and that for 
the FAM-QRPA is given in Ref. [lOj. The implementa- 
tion of a traditional matrix formulation QRPA method 
(MQRPA) used in this work follows that of Ref. 0. In 
this section, we recapitulate the method and define all 
necessary quantities. 

The variation of the total DFT energy defined through 
an energy density £{p, k, k*), with respect to the particle 
and pairing density matrices, p — V*V"^ and k — V*U'^ , 
results in the HFB equations 
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Efi are the quasi-particle energies, (C/^, V^) are the two- 
component HFB quasi-particle vectors, and the chemical 
potential A is introduced to conserve the average particle 
number. 

The QRPA equations for the mode amplitudes (AT^i/, 
y^i/) and excitation energies a; can be written in a matrix 
form as: 
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with matrices A and B coming from second variational 
derivatives of £[p,k,k*] with respect to p and k. In 
MQRPA, Eqs. ^ are solved by means of the explicit 
diagonalization, and the strength function corresponding 
to the one-body operator F is subsequently computed. In 
our results, strength functions calculated with MQRPA 
are smeared with a Lorentzian-averaging function hav- 
ing a width r = 27. Such an averaging can be associated 
with complex QRPA frequencies lo + i"f that are intro- 
duced in the context of stren gth function technique with 
schematic interactions jll. l36l - [38| . In fact, strength func- 
tions obtained in such a way do not require knowledge 
of individual RPA eigenvalues; the summation over the 
RPA spectrum is replaced by integration over energy (see 
also Refs. [ll,[ill). 

Following the earlier applications of FAM [H, [H-d^l , 
we solve the QRPA problem in the presence of a one- 
body external perturbation f of a given frequency w. In 
this case, Eq. ([3]) can be rewritten in an alternative way 
M: 



(E^ + E, - u) X^, + 5Hf,{Lo) = Fll 



(4) 



(5) 



{E^ + E, + Lo) V + 5H^^A^) = F'^l 

where the complex antisymmetric matrices 

5F20(a;) = UHh{uj)V* - VUh{ujYU* 
-V^JK*{uj)V* + WSA{uj)U*, 

dH^^iu) = U'^Sh{ujfV - V'^5h{u)U 
~V'^SA{uj)V + U'^SA{uj)*U, 

are defined in terms of the non-Hermitian variations 

Sh{uj) = {h[pr,, Kr,, Kjj] - h[p, K, K*])/??, 

SA{oj) ={A[p^,K^]-A[p,K])/r,, (6) 

Miiu) =(A[p,„S^]-A[p,K])/77, 

where 77 is a small parameter to numerically expand den- 
sities up to the first order. The non-Hermitian density 
matrices in ([5]) are: 



p,,^ iV + r]U*X*)*{V + r,U*Yy 



* "V^*\t 



Kr, = -([/ + vV*Y) {V + r]U*X*) 
p,,^ (V + ijU*Y)*{V + 'qU*X*f , 
R,,^ -{U + riV*X*){V + rjU*Y)\ 



(7) 



We note that the above density matrices depend on the 
external field F through the QRPA veetors {X, Y). 

The FAM-QRPA equations (|4|) can be formally solved 
with respect to X^^, Y^,^: 
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Since the matrices SH'^'^{lu) and 5H^'^{uj) linearly depend 
on Xfj^i, and Y^ij, a self-consistent iterative scheme needs 
to be adopted to find the QRPA amplitudes. In essence, 
FAM replaces the calculation and diagonalization of the 
large QRPA matrices A, B with a much simpler proce- 
dure of calculating SH^^{uj) and SH^^{uj), and solving 
Eqs. (HI at desired values of w. To guarantee that the 
FAM-QRPA solution carries a finite strength function 
for every value of w, we take a; — > oj + 17 with a small 
imaginary part 7. As we shall see later, such a choice cor- 
responds to a Lorentzian smearing of F = 27, except for 
the vicinity of w = [l| . It is worth noting that the FAM 
implementation of QRPA is straightforward and EDF- 
independent, as the same HFB procedures that define 
the fields h,A in terms of particle and pairing densities 
are also used in FAM-QRPA calculations. 

Ideally, a self-consistent iterative FAM-QRPA proce- 
dure should converge rapidly and the result should not 
depend on ry if its value is small enough. In practice, 
a direct iteration of ([5]) diverges in most cases, espe- 
cially when 7 is small and/or when lo is close to the 
true QRPA root. Indeed, when 7 — )► 0, the left-hand 
side of Eqs. (|3]) becomes singular; hence, instabilities are 
expected around QRPA roots. To guarantee numerical 
stability, one has to resort to a procedure which 'mixes' 
the solutions from previous and next iterations. To this 
end, the conjugate gradient method and its derivatives 
were utilized in coordinate-space applications |26l. l28i - l30l |. 
In this work, based on the HO exp ansion technique, the 
modified Broyden's procedure [SJ, |3g has been adopted 
and turned out to yield stable results while providing 
excellent computational performance. For a system of 
linear equations (J4|), Broyden iterations exhibit the Q- 
superlinear convergence [39[. 

In the Arnoldi (or Lanczos) diagonalization scheme, 
QRPA iterations start from a pivot vector that has its el- 
ements set to the matrix elements of F. This guarantees 
that all odd-power energy-weighed sum rules (EWSR) 
are conserved. The number of states found is equal to the 
number of iterations, and these states (usually different 
from QRPA modes) are used to construct the strength 
function. On the other hand, by solving the linear re- 
sponse equation ^ with a fixed w is that one obtains 
solutions at required energies. In this way the individual 
QRPA states can be further accessed, if needed. 

Results- To benchmark FAM-QRPA for deformed nu- 
clei and check its performance, we carried out FAM and 
MQRPA [U calculations using the SLy4 EDF ^ and 
a contact volume pairing with a 60 MeV cutoff with re- 
spect to the reference single-particle energies [41|. The 



pairing strength was chosen to reproduce the experimen- 
tal neutron pairing gap of ^^"Sn. All HFB calculations 
were performed with the DFT code HFBTHO [4l| that 
solves the Skyrme-HFB equations in the HO basis, as- 
suming axial and reflection symmetries. 

As discussed in Ref. [21| , MQRPA calculations arc sub- 
ject to two truncations. The first truncation pertains 
to the maximum rank of the QRPA Hamiltonian matrix 
that, in our case, should not exceed 20,000. To this end, 
one neglects all canonical states with single-particle en- 
ergies greater than some cutoff value. The second trun- 
cation is made by excluding those QRPA quasiparticle 
pairs that have occupation probabilities less than some 
small critical value v^^i^. or larger than 1 — Wcrit- 

The benchmarking calculations have been carried out 
for the monopolc isoscalar (IS) and modified isovector 
(IV) response operators: 



/^^ = fE^ 



rIV 



7 ^ AT ^ 



(9) 



This choice makes coupling between IS and IV monopole 
vibrations small [42]. At each value w, the imaginary 
part has been set to 7 = 0.5 MeV, and the FAM strength 
function has been calculated according to [30| 

Sif, io) = — ^Im Tr \f{UXV^ + V*Y^U^)\ . (10) 

Here, the external field in Eq. (JH) is given by i^ ~ a/, 
where a is a parameter with dimension [a] = [F][f]^^. 
Since for very small values of 77 the QRPA amplitudes X 
and Y are proportional to a, S{f, to) is independent of a. 
Using a complex frequency uj + ij, the resulting S{f, co) 

Sissesses the crossing symmetry S{f,uj) = —S{f,~Lu) 
; hence S{f,ijj = 0) = is guaranteed. The strength 
function obtained in MQRPA has been computed by av- 
eraging QRPA diagonalization results with F = 27 = 1 
MeV. In the fully self-consistent framework, the MQRPA 
and FAM-QRPA results should be identical. 

The MQRPA and FAM-QRPA J'^ = 0+ strength func- 
tions are compared in Fig. [T] for an oblate minimum in 
^'^Mg. (In the prolate ground state of this nucleus, paring 
correlations vanish.) To include the whole available space 
of canonical wave functions in MQRPA, the results shown 
in Fig. [1] were obtained using a relatively small single- 
particle basis corresponding to A'sh = 5 HO shells. It is 
clearly seen that both methods yield practically identi- 
cal results. Increasing the number of basis states rapidly 
increases the scale of the MQRPA scheme. For example, 
with 20 HO shells and Ucrit = lO^"', QRPA matrices reach 
dimension 32,039, requiring 16.4 GB memory. Lowering 
the canonical cutoff to Udit = lO"^'' results in the matrix 
rank 211,159, or 713.41 GB of memory. In contrast, the 
memory required by FAM-QRPA, using the full space of 
20 HO shells and without any truncations on a QRPA 
level, is a modest 0.572 GB. 

The accuracy of any QRPA implementation can be as- 
sessed by its ability to handle the spurious (zero-energy) 
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FIG. 1: (color online) The isoscalar (blue, dashed line) 
and isovector (red, solid line) nionopole strength function in 
oblate-deformed and paired minimum of Mg obtained in 
MQRPA within the full N^i, = 5 HO space (ucrit=0) and 
FAM-QRPA (circles). The inset shows the presence of 0"*" 
spurious mode at low-energy at three different values of Ucrit 
in MQRPA. 
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modes. In general, QRPA solutions should be properly 
orthogonalized against spurious modes by means of a 
well-established procedure [2a, [23] • For the monopole 
case presented in Fig. [I] the huge configuration space 
of FAM-QRPA seems to be sufficient to remove the 0^ 
spurious modes associated with particle-number break- 
ing almost exactly. This is indeed seen in the inset 
of Fig. [T] which displays the low-energy IV monopole 
strength in MQRPA for several values of Wcrit- The 
MQRPA response corresponding to i)„it=0 (full HO 
space) shows the single low-energy peak carrying the 
strength of about 2-10~^ e^fm^/MeV, and the low-energy 
FAM-QRPA strength is of similar magnitude. 

Figure[2]demonstrates that the FAM results practically 
do not depend on the choice of the parameter rj entering 
the numerical derivatives in Eq. ^ for quite a large 
range of values of 77 from 10~^tolO~^. Usually, the FAM 
solution is reached fairly quickly, within 10-40 iterations 
assuming that the maximum difference between collective 
amplitudes in two consecutive iterations is less that 10~^. 

In order to demonstrate the ability of FAM-QRPA to 
treat heavy, deformed, and superfluid nuclei, Fig.|3]shows 
the monopole strength distributions for the ground state 
and fission isomer in ^*'^Pu obtained with a large con- 
figuration space of TVsh = 20. As it is well known [43j . 
due to its large deformation, the IS monopole strength 
distribution splits into two components. Here, as well as 
in other cases considered in this work, about 99% of the 
EWSR is exhausted when integrating up to a;=50McV. 
Both examples nicely illustrate the applicability of FAM- 
QRPA to the local QRPA approach used in the context 



FIG. 2: (color online) The performance of the FAM-QRPA 
algorithm applied to the case of Fig. [l] for different values of 
77 in the frequency range of 24<aj<32 MeV. Top: number of 
Broyden iterations. Bottom: relative accuracy. 
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FIG. 3: (color onUne) IS (sohd line) and IV (dotted line) 
monopole strength in the deformed ground state of ^''"Pu and 
its superdeformed fission isomer in FAM-QRPA with Vah=20. 



of the large-amplitude collective motion [19| . 

Conclusions- In this work we applied the FAM-QRPA 
method to describe the multipole strength in deformed 
and superfluid nuclei. The calculations have been pre- 
sented for IS and IV monopole modes. We first bench- 
marked FAM-QRPA against the MQRPA approach of 



Ref. 12l| and obtained excellent agreement. As compared 
to the standard diagonalization method, FAM-QRPA of- 
fers excellent performance, both in terms of memory 
and speed. Including all the fields (both time-even and 
time-odd) required by the fully self-consistent QRPA, the 
memory requirement for the FAM-QRPA module built 
on the top of the DFT solver HFBTHO does not ex- 
ceed 572 MB. This enables us to handle axially-deformed 
heavy nuclei without imposing any truncation on the 
QRPA level. The self-consistency of FAM-QRPA, to- 
gether with very large pairing windows employed, results 
in a practical decoupling of the 0"'' spurious modes asso- 
ciated with the particle-number symmetry breaking. 

A new efficient and robust method for the iterative 
solution of FAM-QRPA equations has been proposed. 
The method, based on the Broyden mixing procedure al- 
ready adopted in DFT solvers |35| , offers a Q-superlinear 



convergence of FAM-QRPA equations. The proposed 
FAM implementation allows fast calculations of multi- 
pole strength for all axially deformed nuclei throughout 
the nuclear chart. The algorithm is especially suited 
for multiprocessor tasks since the QRPA strength dis- 
tribution S{uj) converges with no more then 40 iterations 
regardless of w. The implementation of the method to 
higher multipolarity modes is in progress. 
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